Smoothing models

Get data - smooth disturbed and control plots together - AGB growth

library(mgcv); library(data.table); library(gratia); library(marginaleffects); library(ggplot2); library(dplyr); library(tidyr); library(zoo); library(slider); library(here)

#####test Paracou####

obs <- read.csv("G:/My Drive/cesab bioforest/Data/aggregated_data_v9.csv")

setDT(obs)

obs=obs[obs$Site=="Paracou"]


setnames(obs, "Plot", "subplot")

obs$plot=gsub("_.", "",obs$subplot )

obs <- obs[!Year%in%1984:1990]

obs <- dcast(
  obs,
  subplot + plot + Year ~ variable,
  value.var = "value"
)


##### upload data 
load("../clim_by_site.rda")

clim_all=rbindlist(lapply(clim_by_site[["Paracou"]], function(x) x$climate[[1]]),
                   idcol = "plot")

#"Paracou", "Mbaiki", "Corinto", "SUAS", "Lesong", "Ulu Muda", "Sungai Lalang", "Tene", "Jari"

baseline <- clim_all[year >= 1981 & year <= 2010]


clim_stats <- baseline[, .(
  mean_tmax = mean(tmax),
  sd_tmax   = sd(tmax),
  mean_vpd   = mean(vpd),
  sd_vpd     = sd(vpd),
  mean_srad  = mean(srad),
  sd_srad    = sd(srad),
  mean_def  = mean(def),
  sd_def    = sd(def)
), by = month]

clim_all <- merge(clim_all, clim_stats, by = "month", all.x = TRUE)

clim_all[, `:=`(
  z_tmax = (tmax - mean_tmax) / sd_tmax,
  z_vpd  = (vpd  - mean_vpd)  / sd_vpd,
  z_srad = (srad - mean_srad) / sd_srad,
  z_def  = (def  - mean_def)  / sd_def
)]

#climate anomalies follow inventories. For Paracou, mean per year before 1996, two years after

clim_all[, year_bin :=
           ifelse(year <= 1995,
                  year,                     # keep original year
                  1997 + ((year - 1996) %/% 2) * 2   # 2-year bins
           )]


census_anom <- clim_all[,.(
  mean_z_tmax = mean(z_tmax, na.rm = TRUE),
  mean_z_vpd = mean(z_vpd, na.rm = TRUE),
  mean_z_srad = mean(z_srad, na.rm = TRUE),
  mean_z_def = mean(z_def, na.rm = TRUE)
),
by = .(plot, year_bin)]

#AGB values before 1991 are not consistent
census_anom=census_anom[year_bin>=1991] 


merged <- merge(
  obs,
  census_anom,
  by.x = c("Year","plot"), by.y=c("year_bin","plot")
)


control=c(1,6,11,13,14,15)

merged$Treatment <- ifelse(merged$plot %in% control, "control", "disturbed")

merged$Treatment <- factor(merged$Treatment)
merged$plot <- factor(merged$plot)
merged$subplot <- factor(merged$subplot)

GAM k=5, k=8 and bs = “fs”

NAs for agb_growth - plots 13,14 and 15 miss data (?)

merged <- merged %>%
  group_by(subplot) %>%
  mutate(
    fitted_GAM5 = predict(
      gam(agb_growth ~ s(Year, k = 5), data = cur_data()),
      newdata = cur_data()
    )
  ) %>%
  ungroup()


p1=ggplot(merged, aes(x = Year, y = fitted_GAM5, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="GAM-k=5") + theme(legend.position = "none")


plot(merged$agb_growth, merged$fitted_GAM5)

p1.1 <- ggplot(merged, aes(Year, agb_growth - fitted_GAM5)) +
  geom_point() +
  ggtitle("GAM residuals")

p1 + p1.1

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_GAM5)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("GAM5 residuals vs Tmax")

ggplot(merged, aes(x = Year, y = fitted_GAM5, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~ subplot, scales = "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="GAM-k=5") + theme(legend.position = "none")

GAM k=8

merged <- merged %>%
  group_by(subplot) %>%
  mutate(
    fitted_GAM8 = predict(
      gam(agb_growth ~ s(Year, k = 8), data = cur_data()),
      newdata = cur_data()
    )
  )%>%
  ungroup()

p11=ggplot(merged, aes(x = Year, y = fitted_GAM8, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~ subplot, scales = "free_y")+
  
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="GAM-k=8") + theme(legend.position = "none")

p11

plot(merged$agb_growth, merged$fitted_GAM8)

p11.1 <- ggplot(merged, aes(Year, agb_growth - fitted_GAM8)) +
  geom_point() +
  ggtitle("GAM residuals")

 
p11.1

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_GAM8)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("GAM8 residuals vs Tmax")

GAM5 bs = “fs”

#test GAM that is done on the full models (bs = "fs")

merged <- merged %>%
  mutate(
    fitted_GAM5_fs = predict(
      gam(agb_growth ~ s(Year, subplot, bs = "fs", k = 5), data = cur_data()),
      newdata = cur_data()
    )
  )

p111=ggplot(merged, aes(x = Year, y = fitted_GAM5_fs, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  theme_bw() +
  facet_wrap(~ subplot, scales = "free_y")+
  labs(x = "Year", 
       y = "agb_growth", title="GAM-k=5 fs") + theme(legend.position = "none")

p111 

plot(merged$agb_growth, merged$fitted_GAM5_fs)

p111.1 <- ggplot(merged, aes(Year, agb_growth - fitted_GAM5_fs)) +
  geom_point() +
  ggtitle("GAM residuals")

p111.1

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_GAM5_fs)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("GAM5 bs = “fs” residuals vs Tmax")

#GAM captures early and late patterns faithfully, best edges behaviour. more variance on early years not linked to the model

GAM8 bs = “fs”

#test GAM that is done on the full models (bs = "fs")

merged <- merged %>%
  mutate(
    fitted_GAM8_fs = predict(
      gam(agb_growth ~ s(Year, subplot, bs = "fs", k = 8), data = cur_data()),
      newdata = cur_data()
    )
  )

p111=ggplot(merged, aes(x = Year, y = fitted_GAM8_fs, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  theme_bw() +
  facet_wrap(~ subplot, scales = "free_y")+
  labs(x = "Year", 
       y = "agb_growth", title="GAM-k=8 fs") + theme(legend.position = "none")

p111 

plot(merged$agb_growth, merged$fitted_GAM8_fs)

p111.1 <- ggplot(merged, aes(Year, agb_growth - fitted_GAM8_fs)) +
  geom_point() +
  ggtitle("GAM residuals")

p111.1

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_GAM8_fs)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("GAM8 bs = “fs” residuals vs Tmax")

#GAM captures early and late patterns faithfully, best edges behaviour. more variance on early years not linked to the model

slider

complete=TRUE (reject incomplete windows, NA at the beggining and the end), window per observation, not true “moving time windows”, window is just row-based.

merged <- merged %>%
  group_by(subplot) %>%
  arrange(Year) %>%
  mutate(fitted_slide_complete = slide_dbl(agb_growth, mean, .before = 2, .after = 2, .complete=TRUE))%>%
  ungroup()

p3=ggplot(merged, aes(x = Year, y = fitted_slide_complete, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales = "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="slide - .before = 2, .after = 2, window per observation")+ theme(legend.position = "none")

p3 

plot(merged$agb_growth, merged$fitted_slide_complete)

p3.1 <- ggplot(merged, aes(Year, agb_growth - fitted_slide_complete)) +
  geom_point() +
  ggtitle("slide residuals")

p3.1

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_slide_complete)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("slide_complete residuals vs Tmax")

complete=TRUE (reject incomplete windows, increased window NA at the beggining and the end), window per observation, not true “moving time windows”, window is just row-based.

merged <- merged %>%
  group_by(subplot) %>%
  arrange(Year) %>%
  mutate(fitted_slide_complete3 = slide_dbl(agb_growth, mean, .before = 3, .after = 3, .complete=TRUE))%>%
  ungroup()

p3=ggplot(merged, aes(x = Year, y = fitted_slide_complete3, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales = "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="slide - .before = 3, .after = 3, window per observation")+ theme(legend.position = "none")

p3 

plot(merged$agb_growth, merged$fitted_slide_complete3)

p3.1 <- ggplot(merged, aes(Year, agb_growth - fitted_slide_complete3)) +
  geom_point() +
  ggtitle("slide residuals")

p3.1

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_slide_complete3)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("slide_complete3 residuals vs Tmax")

**complete=TRUE, window per Year (NAs can be in the middle - anywhere Year gaps break the window)

merged <- merged %>%
  group_by(subplot) %>%
  arrange(Year) %>%
  mutate(fitted_slide_completeYear = slide_index_dbl(agb_growth, Year, mean, .before = 2, .after = 2, .complete=TRUE))%>%
  ungroup()

p33=ggplot(merged, aes(x = Year, y = fitted_slide_completeYear, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales= "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="slide - .before = 2, .after = 2, window per year")+ theme(legend.position = "none")

p33 

plot(merged$agb_growth, merged$fitted_slide_completeYear)

p33.1 <- ggplot(merged, aes(Year, agb_growth - fitted_slide_completeYear)) +
  geom_point() +
  ggtitle("slide residuals")

p33.1

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_slide_completeYear)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("slide_completeYear residuals vs Tmax")

**complete=TRUE, window per Year - increase window (more NAs)

merged <- merged %>%
  group_by(subplot) %>%
  arrange(Year) %>%
  mutate(fitted_slide_completeYear3 = slide_index_dbl(agb_growth, Year, mean, .before = 3, .after = 3, .complete=TRUE))%>%
  ungroup()

p33=ggplot(merged, aes(x = Year, y = fitted_slide_completeYear3, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales= "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="slide - .before = 3, .after = 3, window per year")+ theme(legend.position = "none")

p33 

plot(merged$agb_growth, merged$fitted_slide_completeYear3)

p33.1 <- ggplot(merged, aes(Year, agb_growth - fitted_slide_completeYear3)) +
  geom_point() +
  ggtitle("slide residuals")

p33.1

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_slide_completeYear3)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("slide_completeYear3 residuals vs Tmax")

complete=FALSE, no NAs, use less points

merged <- merged %>%
  group_by(subplot) %>%
  arrange(Year) %>%
  mutate(fitted_slide = slide_dbl(agb_growth, mean, .before = 2, .after = 2, .complete=FALSE))%>%
  ungroup()

p333=ggplot(merged, aes(x = Year, y = fitted_slide, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales = "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="slide - .before = 2, .after = 2 - partial windows at edges")+ theme(legend.position = "none")

p333 

plot(merged$agb_growth, merged$fitted_slide)

p333.1 <- ggplot(merged, aes(Year, agb_growth - fitted_slide)) +
  geom_point() +
  ggtitle("slide residuals")


p333.1

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_slide)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("slide residuals vs Tmax")

complete=FALSE, increased window no NAs, use less points

merged <- merged %>%
  group_by(subplot) %>%
  arrange(Year) %>%
  mutate(fitted_slide3 = slide_dbl(agb_growth, mean, .before = 3, .after = 3, .complete=FALSE))%>%
  ungroup()

p333=ggplot(merged, aes(x = Year, y = fitted_slide3, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales = "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="slide - .before = 3, .after = 3 - partial windows at edges")+ theme(legend.position = "none")

p333 

plot(merged$agb_growth, merged$fitted_slide3)

p333.1 <- ggplot(merged, aes(Year, agb_growth - fitted_slide3)) +
  geom_point() +
  ggtitle("slide residuals")


p333.1

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_slide3)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("slide3 residuals vs Tmax")

polynome

use ALL years to fit the curve at each year

merged <- merged %>%
  group_by(subplot) %>%
  mutate(fitted_poly2 = fitted(lm(agb_growth ~ poly(Year, 2), data = cur_data(),na.action = na.exclude)))%>%
  ungroup()

p4=ggplot(merged, aes(x = Year, y = fitted_poly2, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales="free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="Quad poly")+ theme(legend.position = "none")

p4

plot(merged$agb_growth, merged$fitted_poly2)

p4.1 <- ggplot(merged, aes(Year, agb_growth - fitted_poly2)) +
  geom_point() +
  ggtitle("Poly residuals")

p4.1

ggplot(merged, aes(mean_z_tmax, agb_growth - fitted_poly2)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("Poly residuals vs Tmax")

root mean square error

rmse_results <- merged %>%
  group_by(subplot) %>%
  summarise(
    RMSE_GAM8   = sqrt(mean((agb_growth - fitted_GAM8)^2, na.rm = TRUE)),
    RMSE_GAM8_fs   = sqrt(mean((agb_growth - fitted_GAM8_fs)^2, na.rm = TRUE)),
    RMSE_GAM5_fs   = sqrt(mean((agb_growth - fitted_GAM5_fs)^2, na.rm = TRUE)),
    RMSE_GAM5   = sqrt(mean((agb_growth - fitted_GAM5)^2, na.rm = TRUE)),
    RMSE_SLIDE  = sqrt(mean((agb_growth - fitted_slide)^2, na.rm = TRUE)),
    RMSE_SLIDE_COMPLETE  = sqrt(mean((agb_growth - fitted_slide_complete)^2, na.rm = TRUE)),
    RMSE_SLIDE_COMPLETE_YEAR  = sqrt(mean((agb_growth - fitted_slide_completeYear)^2, na.rm = TRUE)),
    RMSE_SLIDE3  = sqrt(mean((agb_growth - fitted_slide3)^2, na.rm = TRUE)),
    RMSE_SLIDE_COMPLETE3  = sqrt(mean((agb_growth - fitted_slide_complete3)^2, na.rm = TRUE)),
    RMSE_SLIDE_COMPLETE_YEAR3  = sqrt(mean((agb_growth - fitted_slide_completeYear3)^2, na.rm = TRUE)),
    RMSE_PolyQuad  = sqrt(mean((agb_growth - fitted_poly2)^2, na.rm = TRUE))
  ) %>% 
  ungroup()

rmse_results <- rmse_results %>%
  mutate(
    best_method = colnames(.)[
      apply(select(., -subplot), 1, function(x) which.min(replace(x, is.na(x), Inf))) + 1
    ]
  )

as.data.table(rmse_results)

    subplot RMSE_GAM8 RMSE_GAM8_fs RMSE_GAM5_fs RMSE_GAM5 RMSE_SLIDE
     <fctr>     <num>        <num>        <num>     <num>      <num>
 1:     1_1 0.6602780    0.6579413    0.6572358 0.6597762  0.6862318
 2:     1_2 0.5440815    0.5416930    0.5386886 0.5440815  0.6040639
 3:     1_3 0.6609940    0.6557081    0.6580231 0.6609940  0.7424553
 4:     1_4 0.5800805    0.5803164    0.5839456 0.5800805  0.5960658
 5:    10_1 0.7849051    0.8456235    0.8455720 0.7825946  0.7973029
 6:    10_2 0.6821700    0.6956048    0.7021774 0.6818336  0.6831134
 7:    10_3 0.7672751    0.7549499    0.7543611 0.7682380  0.7282196
 8:    10_4 0.7838802    0.8222049    0.8209406 0.7819798  0.7669645
 9:    11_1 0.7357129    0.7247332    0.7254122 0.7351039  0.7745976
10:    11_2 0.7142709    0.6952320    0.6938137 0.7142709  0.6943462
11:    11_3 0.8031754    0.7794432    0.7858443 0.8031754  0.8220464
12:    11_4 0.8146776    0.7937353    0.7936181 0.8146776  0.7326208
13:    12_1 0.8704537    0.8754618    0.8790203 0.8701286  0.7861002
14:    12_2 0.7417773    0.7371563    0.7388891 0.7417773  0.7541769
15:    12_3 0.9156824    0.9188433    0.9284907 0.9150931  0.8331223
16:    12_4 0.8799091    0.8639448    0.8651460 0.8799091  0.8558761
17:    13_1 0.9668969    0.9638537    0.9641385 0.9692105  0.8783367
18:    13_2 0.9318155    0.9050310    0.9069623 0.9318155  0.8934179
19:    13_3 0.6107574    0.7332620    0.7313225 0.6181148  0.6323810
20:    13_4 0.7650064    0.7648394    0.7757871 0.7651549  0.6167025
21:    14_1 0.7517864    0.7661681    0.7688123 0.7533733  0.7121477
22:    14_2 0.4812247    0.5133870    0.5209513 0.4831050  0.4409424
23:    14_3 0.6529933    0.6871827    0.7006244 0.6535501  0.5654100
24:    14_4 0.6397579    0.6414662    0.6594459 0.6406264  0.5366713
25:    15_1 0.7796617    0.7620071    0.7948504 0.7796617  0.5999784
26:    15_2 0.6667609    0.6523123    0.6640424 0.6667603  0.6222884
27:    15_3 0.6817246    0.6731469    0.6886229 0.6823092  0.7056651
28:    15_4 1.0042997    0.9972733    1.0010379 1.0056790  0.7683992
29:     2_1 0.3835665    0.6019324    0.6088594 0.6227459  0.5444585
30:     2_2 0.4542383    0.5151127    0.5189413 0.4572000  0.4443997
31:     2_3 0.6446780    0.7004875    0.7055906 0.6427212  0.6508805
32:     2_4 0.6203826    0.6101165    0.6116113 0.6203826  0.5947166
33:     3_1 0.7336957    0.7332287    0.7369383 0.7336957  0.6931810
34:     3_2 0.4190214    0.4154055    0.4170717 0.4190214  0.4327896
35:     3_3 0.5668004    0.9629471    0.9471480 0.7456780  0.6978298
36:     3_4 0.4789723    0.4931243    0.4969832 0.4785594  0.4325711
37:     4_1 0.3263275    0.3896101    0.4087643 0.3281933  0.3346153
38:     4_2 0.3738629    0.3762309    0.3765367 0.3736529  0.4240364
39:     4_3 0.6534133    0.7128752    0.7418038 0.6552418  0.6276590
40:     4_4 0.5099975    0.5735599    0.5760701 0.5119338  0.4938520
41:     5_1 0.6145998    0.5986149    0.5999166 0.6145998  0.5774668
42:     5_2 0.4819202    0.4776803    0.5087922 0.4826530  0.4603537
43:     5_3 0.6389234    0.6381846    0.6493634 0.6389234  0.7019902
44:     5_4 0.7110561    0.7032076    0.7142058 0.7110561  0.6971347
45:     6_1 1.0523562    1.0247354    1.0271944 1.0525408  1.0084021
46:     6_2 0.7472778    0.7234551    0.7500571 0.7472778  0.7026307
47:     6_3 0.8233057    0.8158325    0.8166728 0.8233057  0.8410777
48:     6_4 0.7529320    0.7841702    0.7913046 0.7293806  0.7669631
49:     7_1 0.6001234    0.5973120    0.5983429 0.6001234  0.6059278
50:     7_2 0.7373637    0.7359204    0.7438676 0.7373637  0.7936859
51:     7_3 0.6102879    0.6101203    0.6237725 0.6102879  0.6375042
52:     7_4 0.6972089    0.6979671    0.7010979 0.6972089  0.6674518
53:     8_1 0.6704648    0.6765673    0.6838572 0.6632378  0.5724463
54:     8_2 0.7272951    0.7864109    0.8147696 0.7216743  0.7005143
55:     8_3 0.6863311    0.6734833    0.6914963 0.6855547  0.6281695
56:     8_4 0.6313004    0.6696157    0.8118365 0.6302679  0.6508777
57:     9_1 0.6122220    0.6115483    0.6141407 0.6122220  0.5848557
58:     9_2 0.9248327    0.9117891    0.9125826 0.9248327  0.9307054
59:     9_3 0.5965190    0.5913531    0.5948557 0.5966956  0.5728574
60:     9_4 0.7001091    0.6947434    0.6963900 0.7001091  0.6837235
    subplot RMSE_GAM8 RMSE_GAM8_fs RMSE_GAM5_fs RMSE_GAM5 RMSE_SLIDE
    RMSE_SLIDE_COMPLETE RMSE_SLIDE_COMPLETE_YEAR RMSE_SLIDE3
                  <num>                    <num>       <num>
 1:           0.6924368                0.6620793   0.6096979
 2:           0.6558556                0.5574081   0.5101693
 3:           0.7302887                0.6020966   0.6554588
 4:           0.6077942                0.5521897   0.5975627
 5:           0.7615287                0.8883787   0.8356390
 6:           0.5688687                0.6474647   0.6968462
 7:           0.6330230                0.7418817   0.7423033
 8:           0.7733975                0.8172911   0.7674571
 9:           0.6619965                0.7246492   0.6749172
10:           0.6526901                0.6654209   0.6619616
11:           0.6402996                0.6816309   0.7486755
12:           0.5623786                0.6018615   0.8335602
13:           0.6936347                0.9333572   0.9233959
14:           0.6736246                0.6983954   0.7707190
15:           0.7329491                1.0766141   0.9516255
16:           0.7499191                0.8705708   0.8769783
17:           0.8533502                0.9836345   1.0034173
18:           0.8904807                0.9827753   0.9184356
19:           0.4058699                0.6067523   0.7570875
20:           0.5866013                0.6967027   0.6962077
21:           0.6788017                0.7558908   0.7013854
22:           0.3884776                0.3321044   0.4360241
23:           0.5026993                0.5559724   0.5021700
24:           0.4515229                0.5204747   0.5531739
25:           0.5157961                0.5710494   0.5901447
26:           0.6637365                0.5446195   0.5488969
27:           0.6263003                0.6668891   0.7158554
28:           0.7643632                0.8526297   0.9105331
29:           0.4984113                0.4239263   0.6383065
30:           0.4632715                0.4124832   0.5162462
31:           0.5823183                0.5416240   0.6876508
32:           0.6498110                0.6627187   0.6278929
33:           0.7493849                0.8052941   0.7191083
34:           0.4018967                0.3380455   0.4263443
35:           0.6317521                0.6565301   0.8870188
36:           0.4409795                0.5441225   0.4908853
37:           0.3718257                0.3777425   0.3163118
38:           0.4424177                0.3714502   0.3782639
39:           0.5018227                0.7070766   0.6928661
40:           0.5429312                0.5322199   0.5430334
41:           0.6168627                0.6321946   0.5896378
42:           0.4977874                0.4988723   0.4748229
43:           0.6431314                0.6411788   0.6142441
44:           0.7800434                0.8189429   0.7006881
45:           0.8832837                0.8802776   1.0247466
46:           0.7221726                0.6071806   0.7111109
47:           0.5719515                0.7388309   0.7834991
48:           0.6417194                0.6385596   0.7701354
49:           0.6295481                0.6167540   0.5876914
50:           0.8030928                0.7717808   0.7557800
51:           0.6110155                0.6849377   0.6332823
52:           0.6693527                0.8425252   0.7439903
53:           0.5828610                0.7311958   0.7286974
54:           0.7502980                0.8807528   0.7414041
55:           0.6672173                0.7487241   0.6689521
56:           0.7189026                0.6477941   0.6703115
57:           0.6333529                0.5900400   0.6687285
58:           0.9564598                0.9876889   0.9497196
59:           0.5418454                0.5428450   0.6260309
60:           0.5723792                0.6445077   0.7093027
    RMSE_SLIDE_COMPLETE RMSE_SLIDE_COMPLETE_YEAR RMSE_SLIDE3
    RMSE_SLIDE_COMPLETE3 RMSE_SLIDE_COMPLETE_YEAR3 RMSE_PolyQuad
                   <num>                     <num>         <num>
 1:            0.6455194                 0.6558982     0.6516454
 2:            0.5862787                 0.5812258     0.5352876
 3:            0.6466898                 0.6380049     0.6561903
 4:            0.6633609                 0.5592650     0.5793231
 5:            0.8333277                 0.9178252     0.7631782
 6:            0.6268988                 0.6319195     0.6665390
 7:            0.6780787                 0.7435297     0.7522101
 8:            0.8092599                 0.8610417     0.8085605
 9:            0.6329370                 0.7101684     0.7102420
10:            0.6053223                 0.6844717     0.6958772
11:            0.6355198                 0.6677569     0.7745037
12:            0.6024907                 0.6721818     0.7902398
13:            0.9079425                 0.9447773     0.8510292
14:            0.7208718                 0.7131542     0.7400901
15:            0.8838141                 1.0542561     0.8982756
16:            0.8477452                 0.8961059     0.8643039
17:            1.0735848                 1.0223583     0.9410003
18:            0.9887751                 1.0009637     0.9016847
19:            0.5391999                 0.5177383     0.6737868
20:            0.6883399                 0.7312721     0.7444854
21:            0.7135058                 0.7992182     0.7276886
22:            0.3785188                 0.3915866     0.4667363
23:            0.4430737                 0.4789432     0.6330926
24:            0.4264619                 0.4759118     0.6219419
25:            0.5141824                 0.5572178     0.7528862
26:            0.5818494                 0.4269522     0.6430975
27:            0.6694601                 0.7084989     0.6579428
28:            0.9534533                 0.9516403     0.9706199
29:            0.6117028                 0.3734045     0.6122367
30:            0.5402852                 0.4642630     0.5254208
31:            0.5735337                 0.4322319     0.6288648
32:            0.6241531                 0.6399917     0.6044826
33:            0.6653395                 0.8296758     0.7236745
34:            0.4007910                 0.3287068     0.4172664
35:            0.7112646                 0.7278038     0.9045353
36:            0.4845344                 0.5325212     0.4836408
37:            0.3615105                 0.3828170     0.3868553
38:            0.3757339                 0.3877203     0.3632610
39:            0.5580313                 0.6158825     0.6354674
40:            0.5669344                 0.5791783     0.5234642
41:            0.6422167                 0.6817175     0.5946510
42:            0.5187962                 0.5126996     0.4721498
43:            0.5935264                 0.5825218     0.6383874
44:            0.7829082                 0.8632153     0.7015761
45:            0.9239384                 0.9040820     1.0070110
46:            0.7087226                 0.7089890     0.7162201
47:            0.5187326                 0.5422192     0.8130082
48:            0.6610041                 0.6945377     0.7735994
49:            0.6075227                 0.5685193     0.5991195
50:            0.6304469                 0.6414605     0.7367826
51:            0.6414886                 0.7206914     0.6102132
52:            0.7951677                 0.8697104     0.6969681
53:            0.7911909                 0.7798006     0.6814875
54:            0.8288955                 0.9389928     0.7636499
55:            0.6546863                 0.7795544     0.6749501
56:            0.7706704                 0.6871497     0.6291198
57:            0.7161823                 0.6123807     0.6117652
58:            1.0215227                 0.9726187     0.9120191
59:            0.6126634                 0.4970348     0.5772581
60:            0.6288684                 0.5906843     0.6963614
    RMSE_SLIDE_COMPLETE3 RMSE_SLIDE_COMPLETE_YEAR3 RMSE_PolyQuad
                  best_method
                       <char>
 1:               RMSE_SLIDE3
 2:               RMSE_SLIDE3
 3:  RMSE_SLIDE_COMPLETE_YEAR
 4:  RMSE_SLIDE_COMPLETE_YEAR
 5:       RMSE_SLIDE_COMPLETE
 6:       RMSE_SLIDE_COMPLETE
 7:       RMSE_SLIDE_COMPLETE
 8:                RMSE_SLIDE
 9:      RMSE_SLIDE_COMPLETE3
10:      RMSE_SLIDE_COMPLETE3
11:      RMSE_SLIDE_COMPLETE3
12:       RMSE_SLIDE_COMPLETE
13:       RMSE_SLIDE_COMPLETE
14:       RMSE_SLIDE_COMPLETE
15:       RMSE_SLIDE_COMPLETE
16:       RMSE_SLIDE_COMPLETE
17:       RMSE_SLIDE_COMPLETE
18:       RMSE_SLIDE_COMPLETE
19:       RMSE_SLIDE_COMPLETE
20:       RMSE_SLIDE_COMPLETE
21:       RMSE_SLIDE_COMPLETE
22:  RMSE_SLIDE_COMPLETE_YEAR
23:      RMSE_SLIDE_COMPLETE3
24:      RMSE_SLIDE_COMPLETE3
25:      RMSE_SLIDE_COMPLETE3
26: RMSE_SLIDE_COMPLETE_YEAR3
27:       RMSE_SLIDE_COMPLETE
28:       RMSE_SLIDE_COMPLETE
29: RMSE_SLIDE_COMPLETE_YEAR3
30:  RMSE_SLIDE_COMPLETE_YEAR
31: RMSE_SLIDE_COMPLETE_YEAR3
32:                RMSE_SLIDE
33:      RMSE_SLIDE_COMPLETE3
34: RMSE_SLIDE_COMPLETE_YEAR3
35:                 RMSE_GAM8
36:                RMSE_SLIDE
37:               RMSE_SLIDE3
38:             RMSE_PolyQuad
39:       RMSE_SLIDE_COMPLETE
40:                RMSE_SLIDE
41:                RMSE_SLIDE
42:                RMSE_SLIDE
43: RMSE_SLIDE_COMPLETE_YEAR3
44:                RMSE_SLIDE
45:  RMSE_SLIDE_COMPLETE_YEAR
46:  RMSE_SLIDE_COMPLETE_YEAR
47:      RMSE_SLIDE_COMPLETE3
48:  RMSE_SLIDE_COMPLETE_YEAR
49: RMSE_SLIDE_COMPLETE_YEAR3
50:      RMSE_SLIDE_COMPLETE3
51:              RMSE_GAM8_fs
52:                RMSE_SLIDE
53:                RMSE_SLIDE
54:                RMSE_SLIDE
55:                RMSE_SLIDE
56:             RMSE_PolyQuad
57:                RMSE_SLIDE
58:              RMSE_GAM8_fs
59: RMSE_SLIDE_COMPLETE_YEAR3
60:       RMSE_SLIDE_COMPLETE
                  best_method

lm

vars <- grep("^fitted_", names(merged), value = TRUE)

results <- data.frame(
  method = vars,
  intercept = NA,
  slope = NA,
  p_value = NA
)



for (i in seq_along(vars)) {
  
  v <- vars[i]
  model <- lm(as.formula(paste0(v, " - agb_growth ~ Year")), data = merged)
  s <- summary(model)
  
  results$intercept[i] <- coef(s)[1, 1]
  results$slope[i]     <- coef(s)[2, 1]
  results$p_value[i]   <- coef(s)[2, 4]
}

results

                       method     intercept         slope    p_value
1                 fitted_GAM5  2.926468e-11 -1.459840e-14 1.00000000
2                 fitted_GAM8  3.746029e-11 -1.375740e-14 1.00000000
3              fitted_GAM5_fs -4.711339e+00  2.350138e-03 0.31111357
4              fitted_GAM8_fs -7.298793e-01  3.640826e-04 0.87387900
5       fitted_slide_complete  8.845040e+00 -4.412354e-03 0.13032049
6      fitted_slide_complete3  1.964098e+01 -9.806315e-03 0.01022311
7   fitted_slide_completeYear -5.540808e+00  2.780665e-03 0.32114410
8  fitted_slide_completeYear3  1.305033e+01 -6.503793e-03 0.05966999
9                fitted_slide  4.140980e+00 -2.064655e-03 0.35270548
10              fitted_slide3  1.826779e+00 -9.074865e-04 0.69530117
11               fitted_poly2 -4.420647e-12  2.205123e-15 1.00000000

temporal signal on residuals

methods <- c(
  "fitted_GAM5",
  "fitted_GAM8",
  "fitted_GAM8_fs",
  "fitted_GAM5_fs",
  "fitted_slide",
  "fitted_slide_completeYear",
  "fitted_slide_complete",
  "fitted_slide3",
  "fitted_slide_completeYear3",
  "fitted_slide_complete3",
  "fitted_poly2"
)

for (m in methods) {
  merged[[paste0("resid_", m)]] <- merged$agb_growth - merged[[m]]
}

pairwise_corr <- function(df, col) {
  wide <- tidyr::pivot_wider(df[, c("subplot","Year", col)],
                             names_from = subplot, values_from = all_of(col))
  M <- as.matrix(wide[, -1])          # remove Year
  C <- cor(M, use = "pairwise.complete.obs")
  mean(C[upper.tri(C)], na.rm = TRUE)
}

temporal_signal <- data.frame(
  method = methods,
  mean_pairwise_corr = sapply(paste0("resid_", methods),
                              function(x) pairwise_corr(merged, x))
)

temporal_signal

                                                     method mean_pairwise_corr
resid_fitted_GAM5                               fitted_GAM5          0.3102782
resid_fitted_GAM8                               fitted_GAM8          0.3110250
resid_fitted_GAM8_fs                         fitted_GAM8_fs          0.3041428
resid_fitted_GAM5_fs                         fitted_GAM5_fs          0.2994613
resid_fitted_slide                             fitted_slide          0.3169307
resid_fitted_slide_completeYear   fitted_slide_completeYear          0.2853856
resid_fitted_slide_complete           fitted_slide_complete          0.2960472
resid_fitted_slide3                           fitted_slide3          0.2972938
resid_fitted_slide_completeYear3 fitted_slide_completeYear3          0.2937522
resid_fitted_slide_complete3         fitted_slide_complete3          0.3368332
resid_fitted_poly2                             fitted_poly2          0.3120015